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In this paper we present the results of numerical studies of the JIMWLK and BK equations with a particular 



emphasis on the universal scaling properties and phase space structure involved. The results are valid for near zero 



impact parameter in DIS. We demonstrate IR safety due to the occurrence of a rapidity dependent saturation 
scale Qs(t). Within the set of initial conditions chosen both JIMWLK and BK equations show remarkable 
, agreement. We point out the crucial importance of running coupling corrections to obtain consistency in the UV. 

Despite the scale breaking induced by the running coupling we find that evolution drives correlators towards an 
asymptotic form with near scaling properties. We discuss asymptotic features of the evolution, such as the r- 
^ ' and A-dependence of Qs away from the initial condition. 

o 

cn 
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2 ! 1 Introduction 

pi • With the advent of modern colliders, from the Tevatron, HERA to RHIC and planned experiments like LHC 
and EIC, the high energy asymptotics of QCD has gained new prominence and importance. The question 
if QCD effects there show universal features that might be calculated without recourse to intrinsically non- 
pH I perturbative methods capable of also dealing with large as has gained new importance. In situations in 
which the experimental probes see large gluon densities at small x the answer to this question has, maybe 
surprisingly, turned out to be positive. The progress has come from a method that resums the "density" 
■ induced corrections exactly while treating other effects, such as the x dependence as perturbative corrections. 
^ I The justification for this is that the density effects induce a generally x dependent correlation length 1/Qs{x) 
which also sets the scale for the coupling, thus giving us a small expansion parameter. There arc several 
properties of the interactions at small x that allow us to perform such a resummation. To explain them let us 
consider the simplest example of a j*A collision, deep inelastic scatting (DIS) of virtual photons off nuclear 
targets, be it protons or larger nuclei. Generalizations to diffraction, pA and AA scattering are possible and 
share many of the key features. 

Thus we are considering a collision of a virtual photon that imparts a (large) space-like momentum q (Q^ := 
—q^ is large) on a nucleus of momentum p. At large energies the only other relevant invariant is Bjorken x 
with x := A leading twist analysis of DIS employing an OPE based on the largeness of allows us 

to express the the cross section entirely in terms of two point functions, the quark and gluon distributions. 
At the same time it becomes possible to calculate their dependence to leading logarithmic accuracy 
(summing corrections of the form (a^lnQ^)") by solving the DGLAP equations. This amounts to summing 
the well-known QCD ladder diagrams. This procedure is self consistent and perturbatively under control 
as long as we stay within the large region. Starting from a large enough and going to larger values, 
we find that the objects we count with the quark and gluon distributions increase in number but they stay 
dilute: their "sizes" follow the transverse resolution scale 1/Q. At small x, however, the growth of the gluon 
distribution is particularly pronounced and any attempt to rcsum corrections of the form (a^ ln(l/.T))") to 
track the x-dependence of the cross sections immediately is faced with extreme growth of gluon distributions: 
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moving towards small x at fixed increases the number of gluons of fixed "size" 1/Q, so that the objects 
resolved will necessarily start to overlap. This is depicted in Fig. ^ At this point we have long since left the 




Figure 1: Gluon distributions grow towards small x. To the left: qualitative growth as from typical DGLAP fits of 
DIS at HERA (eA) for different Q^. To the right: Growth of gluon (parton) numbers with and ln(l/s). Going to 
sufficiently small x at fixed partons start to overlap. 

region in which an OPE treatment at leading twist is meaningful. Instead of single particle properties like 
distribution functions we need to take into account all higher twist effects that are related to the high density 
situation in which the gluons in our target overlap. At small x we are thus faced with a situation in which 
we clearly need to go beyond the standard tools of perturbative QCD usually based on a twist expansion. 

There are different ways to isolate the leading contributions to the cross section at small a; in a background of 

l2Hr22 23 ,24j. We will use a physically quite intuitive picture first used in the formulation of the McLerran- 
Venugopalan model, a formulation that was originally given in the infinite momentum frame of the nuclear 
target. Probing this target at small x implies a large rapidity difference r = ln(l/a;) between projectile and 
target, corresponding to a large relative boost factor e'^. Knowing from the above considerations that the 
interaction will be dominated by gluonic configurations we try to isolate the boost enhanced ones. Looking at 
the gluon field strength of the target for that purpose we find that components are boost enhanced, while 
all others arc left small and can, for our purposes, be set to zero. At the same time we find a strong Lorentz 
contraction, which wc take to be in x~ direction, and a time dilation, correspondingly in .t+ direction. With 
such a field strength tensor there we are left with only one important degree of freedom which, by choice of 
gauge, can be taken to be the + component of a gauge field. Taking into account time dilation and Lorentz 
contraction the gauge field can then be written as 

A^b + SA with 6'^-=0, b+ ^ f3{x)6{x-) . (1) 

where the leading contribution 6+ is independent and Lorentz contracted to a (5-function. Both of these 
are to be taken to be true with a resolution imposed by the rapidity separation r = ln(l/a;). Note that 
mathematically we can always trade a gauge field that has only a single component for a path ordered 
exponential along the direction that picks up this component: 

b+=i{d+U)U^ U =Pexp--i Jdz-b+{z-,x,0) (2) 

If multiple cikonal interactions are relevant as the high energy nature of the process would suggest, we expect 
these path ordered exponentials to be the natural degrees of freedom. 

Ignoring the small fluctuations SA for a moment we have a picture in which the 7* cross section arises from 
a diagram in which the photon splits into a qq pair which then interacts with the background field. Due to 
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the S like support of 6+ the qq are not deflected in the transverse direction during that interaction. This just 
reflects the largeness of the longitudinal momentum component of these partons at small x. This is shown 
in Fig. 2(a) To be sure, the physics content is not frame dependent although it is encoded differently in. 



say the rest frame of the target. There, wc see neither Lorentz contraction nor time dilation. However the 
scale relations are preserved: the photon splits into a qq pair far outside the target and its p~ is so large that 
typical variations of the target are negligible during the interaction. As a consequence the probe is not 
deflected in the transverse direction, picking up any multiple interactions with (gluonic) scattering centers 



as it punches though the target. This is shown in Fig. 2(b) 



Lorentz contracted to ''d[x 



: huge 



' .41/3; finite 





(a) 



(b) 



Figure 2: qq pairs interacting with a nuclear target at small x. (a) shows the situation in the targets infinite 
momentum frame with fully Lorentz-contracted target fields, (b) shows the situation in the target rest frame in which 
the j*qq vertex is far outside the target at a distance proportional to 1/x — e"^ . 



That multiple interactions are of relevance immediately becomes obvious once we try to calculate such 
diagrams with a background field method. To this end we calculate the diagram shown in Fig. [5] in the 
background of a field of the type shown in Eq. ignoring the small fluctuations. This immediately leads 
to the expression 

amsixBi,Q') = Im ™<^Q^™. = Jd\ \^\r'Q') J <fh {I^^—^^) (3) 

where r = x — y corresponds to the transverse size of the qq dipolc and b ~ (x + y)/2 to its impact parameter 
relative to the target. As shown, the leading small x contribution factorizes into a wave function part (the 
factor IV-'^K^^Q^)) a-iid the dipolc cross section part 

f ^ tr(l - U^Ul) 

adipolc(T = ln(l/.T),r2) J <fb ^)r=ln(l/.) (4) 

The former contains the "f*qq vertices and can, alternatively to the background field method, be calculated 
in QED. They turn out to be linear combinations of Bessel functions A'o,i corresponding to the polarizations 
of the virtual photon. It contains all of the direct Q^-dependence of the cross section. 

As anticipated above, the dipole cross section part carries all the interaction with the background gluon 
fields in terms of the path ordered exponentials. The averaging (. . .)x is understood as an average over the 
dominant configurations at a given x. It contains all the information about the QCD action and the target 
wave functions that are relevant at small x. 

It is fairly clear from the above discussion, that the isolation of the leading 6+ contribution that defines this 
average is resolution dependent idea: as we lower x additional modes, up to now contained in 6 A of Eq. ^ 
will take on the features of 6"*". They will Lorentz contract and their x"^-dependence will freeze. Accordingly 
the averaging procedure will have to change. If we write the average as 

(. ..)r = j D[b+] . . . Wr[b+] or, equivalcntly (. . .)^ = j D[U] . . . Zr[U] , (5) 
the weights W^t-[6+] or Zr[U] will, by necessity be r = ln(l/a;) dependent. 
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Indeed, while the weights themselves do contain nonperturbative information, the change of these weights 
with T can be calculated to leading log accm'acy, i.e. to accuracy as ln(l/a;). Assuming we know, say Zt[U] 
at a given tq, this is done in a Wilson renormalization group manner by integrating over the fluctuations SA 
around b'^ between the old and new cutoffs tq and ri. Taking the limit St = ti — tq wc then get a 
renormalization group equation for Zt[U]. 

Because we integrate over 6 A in the background of arbitrary b of the form Eq. the equation treats 
the background field exactly to all orders and captures the all nonlinear effects also in its interaction with 
the target. The resulting RG equation is functional equation for Zr[U] that is nonlinear in U. It sums 
corrections the the leading diagram shown in Eq. (O in which additional gluons are radiated off the initial 
qq pair as shown schematically in Fig. |3| All multiple eikonal scatterings inside the target (the shaded areas) 
are accounted for. 

target 
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exponentially magnified 
w.r.t. physical Icngtli scales 




Figure 3: RG corrections to the average over background field configurations. The shaded areas contain the contri- 
butions subsumed into the averaging procedure at successive values of r = ln(l/a::). 



The result of this procedure is the Jalilian-Marian+Iancu+McLerran+Weigert+Leonidov+Kovner equa- 
tion | 13lll4l[T51ll6lll8lt21„22...23il21) [the order of names was chosen by Al Mueller to give rise to the acronym 
JIMWLK, pronounced "gym walk"]. This equation and its limiting cases will be discussed in detail in the 
next section. 

Let us close here with a few phenomenological expectations for the key ingredient of our phenomenological 
discussion, the dipole cross section. Clearly, we expect the underlying expectation value of the dipole operator. 



xy 



=ln(l/x) 



(6) 



4 



typeset with I^TeX on February 1, 2008, 22:08 



to vanish outside the target, so that the impact parameter integral essentially samples the transverse target 
size and thus scales roughly with A^l'^ . For central collisions, i.e. deep inside the target we expect the 
dependence on the dipole size r = a; — y to interpolate between zero at small r and one at large r. This 
is the idea of color transparency. The simplest model to fit this requirement (which should be good where 
edge effects can be expected to be small, i.e. for large nuclei) would be to parametrize the b-dcpendcnce 
with a ^-function and have the r dependence interpolate between and 1 in a way that the transition is 
characterized by a single scale Qa usually called the saturation scale. 

In such a model, the h integral will only provide a normalization factor. The idea that the scale in the r 
dependence carries at the same time all the r = ln(l/a;)-dcpendcncc according to 

^y = N{{x-yfQa(Tf) (7) 

has actually been used by Golec-Biernat and Wiisthoff |25ll2f)ll77] to provide a new scaling fit of HERA F-2, 
data that is remarkably good. This is particularly true in view of the fact that we are precisely not talking 
about large nuclear targets. A slight generalization of this scaling ansatz has also been used to connect DIS 
data on nuclei and protons based on the simple idea underlying the McLcrran-Vcnugopalan model that Qa 
should increase like 

{QDa = {Q% (8) 

as the additional gluons encountered in additional nucleons on a straight line path through a larger target 
add in additional color charges that decorrelate eikonal scatterers on shorter and shorter distances. Despite 
the fact that the data do not really cover x- and Q^-ranges where such scaling arguments should reliably 
hold, the agreement is surprisingly good.^ It would appear that the data available show a remarkable degree 
of scaling even in cases that are at best at the edge of regions where we would expect the underlying ideas 
to fully apply. 

In the next section will discuss evolution equation and their relation to scaling features as discussed above. 
As the ideas there were relying on central collisions we will find that the evolution equations lead to such 
features there. Indeed they turn out not to be valid in very peripheral regions. 



2 A short review of the evolution equations 

After the qualitative introduction in the previous section, we will simply state the JIMWLK equation and 
then discuss its key features and limiting cases in order to connect its properties with the phcnomcnological 
discussion above. The JIMWLK equation reads 

drZ,\V\^~y^lxty^%Z,\V\ (9) 

and was first cast in this compact form in [21]. Again, are the Wilson line variables describing the 
kinematically enhanced degrees of freedom, are functional version of the left invariant vector fields (c.f. 
Lie derivatives) acting on Uy according to i^^Uy = Uxt^S^xy (for more details see |221) and 



Xxy ■ — ^2 ^ 



[The indicates the adjoint representation.] Zt[U] is the "weight functional" that determines correlators 
0\U] of U fields according to 



{0[U]), J D[U]0[U]Zr[U] (11) 
D[U] in Eq. Hll|) is a functional Haar- measure. 

^This fit was attempted by the authors of the present paper together with Andreas Preund and Andreas Schafer with the 
idea of checking how far presently available data failed to scale — the results were much closer to scaling than expected 
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Eq. (jSJ itself is of the form of a (generalized) functional Fokker-Planck equation. In particular the operator 
on its right hand side, 

H..:=\iyZx^J'y^yi (12) 

is indeed self adjoint and positive definite and may therefore be called the Fokker-Planck Hamiltonian of the 
small X evolution. 

As a functional equation, Eq. Q is equivalent to an infinite set of equations for ?i-point correlators of U 
and U'' fields in any representation of SU{Nc). These will form a coupled hierarchy due to the [/-dependent 
nature of i/pp- The evolution equation for a given correlator of [/-fields is easily obtained by multiplying 
both sides of Eq. Q and integrating the result with a functional Haar measure. Using the self-adjoint nature 
of Hpp and the fact that the Lie derivatives do not interfere with the Haar measure we immediately obtain 



dr{0[U])r = {{-T:^'^ZxZ''y^'^lO[U]))r . (13) 



The only step remaining in the derivation of the evolution equation for 0[U] is to evaluate the Lie derivatives 
as they act on 0[[/]. The nonlinear nature of H^p implies that the expectation value on the r.h.s. will contain 
more U fields than the l.h.s. and thus a new type of correlator of U fields. To determine {0[U])r we thus 
will need also the evolution equation of this new operator in which the same mechanism will couple in still 
further U fields. Continuing the process we end up with an infinite hierarchy of equations. 

As an example for 0[U] with immediate phenomenological relevance we consider the two point function of 
the dipole operator 

tr(l -UlUy) 



Nr,.y := {N^y)r N^y := ^ " (14) 



[Below we will often suppress the explicit r dependence.] Eq. immediately leads to 

a N f (x — v)'^ - ^ ^ ^ ^ 

dr{N.y)r = ^ j'^' \x-zY{z-yY ^^'^'^ + ~ ^^^^ 

Clearly the r.h.s. depends on a 3-point function containing a total of up to 4 [/^t' factors, that in general 
does not factorize into a product of two 2-point correlators (or linear combinations thereof): 

{N^.N,y)^^{N^,)^{N,y)^. (16) 

To completely specify the evolution of {Nxy)T we therefore also need to know (^Nxz Nzy)^- The latter, in its 
evolution equation, will couple to yet higher n-point functions and thus we are faced with one of the infinite 
hierarchies of evolution equations that are completely encoded in the single functional equation . 

Because of the fact that Eq. (jS)) encodes information on the evolution of all possible n-point functions any 
generic statement that can be derived without reference to a particular correlator is of particular importance. 
Interestingly enough, and despite its functional nature it allows us to deduce already some of the main features 
of the evolution process. First of all, Eq. (j^ is of Hamiltonian form with a semi-positive-definite Hamiltonian 
and describes a diffusion process in which r plays the role of (imaginary) time. From this observation alone, 
we get very strong statements without any additional work: the eigenmodes of this Hamiltonian, Z„.t-[[/], 
possesses positive real eigenvalues A„ and a generic solution will evolve as a superposition of such modes in 
the form 

Zr[U] = ^a„e-^"(^-^°)Z„,,J[/] (17) 

n 

where the sum is a formal sum over all modes in the spectrum. If there is an eigenvalue Aq = 0, this will 
be the only one not damped away exponentially with e"'^"'^ at large r. Indeed Zt^\=o[U] = 1 furnishes such 
a candidate and it can be argued to be unique |22j . This solution features as a fixed point of the evolution 
equation and, since all other contributions naturally die away, it is attractive. This remains true although 
Hpp is scale invariant and the spectrum, labelled by n in the above, is continuous. 

At this fixed point, all correlators will be fully determined by what remains in the averaging procedure after 
Zr[U] 1, namely the Haar measure. This however will average everything to zero unless the transverse 
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locations of the factors in the correlator coincide at least pairwise and for each distinct transverse point 
remaining there is a singlet piece in the decomposition of the U factors into irreducible representations. This 
immediately implies that all correlation lengths at r = oo vanish: evolution must erase them at t grows. 
What does that mean for our standard example the two point or dipole correlator? As we have discussed, 
we expect it to interpolate between zero at small separation and one at large separation, with the transition 
say at a scale l/Qs{xo) at some finite xq. 1/Qs{xq) is the correlation length visible in this correlator. As 
we lower x the correlation length will shrink to zero and we expect the following qualitative trend in its x 
evolution: 




Figure 4: Generic evolution trend for a one scale dipole correlator N^y plotted against |a; — y| as 2: ^ 0: the curves 
move towards the left as r = ln(l/a;) increases, the asymptotic form being a step function at the origin 



There is, in principle, a way to implement Eq. ^ numerically. Eq. (j^ is a Fokker-Plank equation and as 
such can be rewritten as a Langevin equation |22[I28| . The latter can be treated numerically with a lattice 
discretization of transverse coordinate space. We will discuss this at length in Sec. 13 

It is useful to first study the qualitative properties of this equation in a simplified setting, the only approx- 
imation scheme known that preserves the main qualitative features depicted in Figure^ This turns out to 
be a \/Nc expansion. Viewed from this perspective Eq. 116() turns into 

(^x. N^y)^ = {N^.)^{N,y)^ + O(i-) . (18) 

This is in fact general: dropping the corrections naturally factorizes all n-point functions into 2-point functions 
and higher correlators decouple for example from Eq. (|15|l . The remainder turns into 

where the t dependence of the dipole function N is implicitly understood. Eq. H19|) was independently 
derived by Ian Balitsky from his "shock wave formalism" |11[I17| that is completely summarized by Eq. lO 
and by Yuri Kovchegov from Mueller's dipole formalism for the case of large nuclei [50]. Indeed the whole 
infinite coupled hierarchy of equations represented by Eq. ^ is reduced in content to this one equation for 
one single function Nzy, called the Balitsky-Kovchegov- or BK-cquation.^ 

The main feature of this equation is that it remains consistent with an interpretation of 1 — Nxy as a correlator 
of Wilson lines (tr(C/j.C/y))/7Vc in that the nonlinearity keeps N between zero and one. Loosely speaking, this 

^Equations for higher correlators of course still exist. They appear as sums with the BK-equation repeated in all terms 
in a way that follows the factorization of the left hand side. For example dri^Nxz N^y'^^ = {drNscz)N;^y + N^^^ {d-rN^y) = 
(right hand side of 'BK)x2:Nzy + A^iiz (right hand side of 'BK)zy, up to terms of order l/N^. 
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is due to a cancellation of the linear and nonlinear terms in N within Eq. (|19|l . In the low density limit, where 
the N'^ term is negligible we are left with the BFKL equation in a space time setting as first demonstrated 
in jl7| . The solution for the equation without the nonlinear term is indeed of the BFKL form 

^r^-y = / ^'^'^0 {{x - y)2)-H-e(--^°)>^('')(r2)-^— iV.„,.„ (20) 

with arbitrary initial condition A^To.ro ^-^d x{^) the BFKL eigenvalue function. This solution is known to 
grow and spread equally to smaller and larger sizes. The first property limits its range of applicability: as 
soon as N has grown enough as to "switch on" the nonlinearities, the result Eq. H20(l will no longer hold. The 
second shows that the linear equation spreads into the infrared. Where that happens we have do abandon 
it altogether. If the nonlinearity is active at short enough distances, i.e. within the perturbative domain we 
expect two effects: the growth of N should be tamed and the infrared should never affect the evolution. 

Numerical studies j29ll8()| using a momentum space representation of the BK equation with a boundary 
condition localized at a single momentum have shown that this is indeed the case. Moreover the solutions 
have been found to approach a scaling regime in which the coordinate and x dependence are coupled according 
to 

Nr.a=y^N{{x-y)^Q,{x)^). (21a) 

There the saturation scale Qs{x) grows like a power of 1/x. However, the initial conditions and hence the 
dipole correlators in these simulations do not satisfy the "color transparency" boundary conditions referred 
to above. 

That this simply means that scaling behavior is according to Eq. H21a() is quite generic was shown in 
There the authors show that the existence of scaling solutions of the type in Eq. (|21a|) , now with with physics 
inspired spatial "color transparency" boundary conditions of the type 

iV(0) = and iV(oo) = 1 , (21b) 

also leads to power like growth of the saturations scale according to 

Qs{x)^[^yQo (22) 

If it is possible to satisfy all of the requirements in Eqns. H21fl . the value of the constant A is determined 
self-consistently through the equation together with the shape of N.^ That the equation indeed possesses 
such solutions has, up to now only been demonstrated numerically and form part of our results below. From 
our simulations it would appear that that single scale initial conditions that satisfy Eq. Ij21bp typically evolve 
into scaling solutions of the form Eq. H21a|l . 

In principle it is possible to extend the scaling argument of to solutions of the JIMWLK equation although 
there one is talking about an infinite hierarchy of different n-point functions and the claim that this would 
naturally fit into a one scale description would sound even more arbitrary than for the BK case. The reason 
why such solutions are important is that they appear to represent the small x asymptotic behavior for a 
large class of initial conditions. Unlike the qualitative features imposed by the infinite energy fixed point of 
Eq. ((5J), this statement is much more directly amenable to experimental verification: if the initial conditions 
realized in nature put us not too far from this asymptotic behavior, we find a unique and universal feature 
in the a;-dependence of any possible target for which the BK equation yields a reasonable description of the 
zero impact parameter contributions. 

It remains important to understand the conditions under which such behavior can be expected: 

First note, that the scaling behavior Eq. H21|l with growing saturation scale Qg is just a special case of the 
generic behavior shown in Fig. 0] Whether we are in a purely one scale situation or in a setting with more 

^The argument in 1311 does not prove the existence of such solutions, it merely translates scaling ansatz + boundary conditions 
into a particular form of the x dependence of the saturation scale given in Eq. 1221 . Indeed the argument would go though even 
without the nonlinearity. In this case, however, the equation is closely related to the BFKL equation, with BFKL eigenvalues 
and power like cigcnmodes. The generic solution then is given in Eq. 1201 and docs not possess solutions with the properties in 
Eq. jn). 
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than one scale enters the scale that controls the transition of N^-^xy 1 at large r = x — y sets the infrared 
scale in the problem as N will not grow beyond 1. For this to happen we need large densities of gluons. Where 
this is guaranteed at some initial to = ln(l/xo), general fixed point considerations and explicit simulations 
show that the IR scales will only get larger as we continue to evolve towards smaller x. Where the initial 
IR scale is not perturbative the derivation of the evolution equations is not valid, however. This limits our 
description to central impact parameters off large enough nuclei. 

Besides these facts there remains a lot to be learned about solutions to the nonlinear evolution equations 
presented in this section. This is the goal of this paper. 

Obviously, the main advantage of Eq. H19|l over the Langevin representation of Eq. 1^ is its simplicity. 
Although, even here, no analytical solutions are available, the numerical task at hand is much more tractable. 
If we can establish that this equation allows us to identify most of the qualitative features of the full evolution, 
we can use it to explore the main features of small x-evolution an only then resort to a numerical treatment 
of Eq. in the most interesting regions to get quantitatively precise results. It is precisely this role in which 
we are going to utilize the BK-equation. 

Our work consists of two main parts: 

• First, we perform a study of the JIMWLK equation, using its Langevin form implemented using a 
lattice discretization of transverse space. 

This is the first coordinate space treatment of either the JIMWLK or BK equations and we are the 
first to show directly that the solutions approach the fixed point as qualitatively predicted in \22\ . c.f. 
Fig.H 

At the same time we will be able to confirm IR safety of the JIMWLK equation via numerical simula- 
tions. 

Still at finite r we will see that scaling solutions emerge that are characterized by a saturation scale 
Qs{x), a feature already found earlier in the context of the BK-equation |;-i()ir29j . 

The initial conditions chosen in this study are almost free of any iVc- violations. We measure how these 
fare under full JIMWLK evolution and find that the equations appear to neither enhance nor dampen 
the extent of these. 

To complement our study of IR stability we also carefully study the UV or continuum limit. Wc will 
show that diffusion into the UV is not tamed. Recalling the properties of the BFKL equation this does 
not come as a complete surprise. However, we will show that for some quantities like the evolution 
rate A,- dr^n{Qs{T)) this fact makes a reliable continuum extrapolation based on JIMWLK a task 
of enormous computational cost. 

Due to the nearly perfect A^c-factorization properties associated with the initial conditions used we are 
then able to complement our JIMWLK studies with BK simulations implemented with same type of 
lattice UV regulator. This allows for a more thorough study of the continuum limit and reveals even 
closer agreement of JIMWLK and BK behavior for the classes of initial condition chosen. We will show 
that the active phase space in both cases extends to up to 7 orders of magnitude above Qg, the physical 
scale of the problem. In such situations one is typically forced to resum large logarithmic corrections 
to get quantitatively reliable results. This then leads us to the second part: 

• Motivated by results found in connection with the BFKL equation we will identify left out running 
coupling effects as the main culprits for the largeness of UV phase space. We perform simulations with 
a coupling running with the size of the parent dipole as as{l/ {x — y)'^) and establish that this removes 
the extraneous phase space and thus resums large phase space corrections. Only with this resummation 
in place the main contributions to evolution originate at distances at and around l/Qs{x) and the 
coupling indeed is of the order of as{Qsix)'^) as anticipated in the underlying physics picture. 

Aqcd and scale breaking enter our simulation as a consequence of this step. This, at least conceptually, 
breaks the asymptotic, geometric scaling limit of the equations. Numerical evidence nevertheless shows 
that the solutions still converge to an approximate asymptotic form that shows almost perfect geometric 
scaling, in agreement with Nevertheless we will show that a clear modification of the x-dependence 
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of Qs emerges, away from the simple fixed coupling (xo/a;)'''-form. Aoff (suitably defined, see below) picks 
up pronounced t dependence and the numerical values change drastically: we find typically a factor of 
two to three compared to the fixed coupling situation. While the necessity of the inclusion of running 
coupling effects had been seen already in estimates of evolution rates and numerical simulations |29II32| 
as well as phcnomcnological fits of HERA data 123 , size and conceptual importance of these corrections 
has not been discussed from the viewpoint taken here. 



3 JIMWLK and BK equations at fixed coupling 
3.1 Simulating the JIMWLK equation 

To be able to solve the JIMWLK equation directly we need a numerical method to implement its content. 
Such a method offers itself after translating the Fokker-Planck equation to a Langevin equation. Such an 
approach has already been suggested in j22| and a derivation has been given in j2g] where we refer the reader 
for details. For the interested, we provide a short exposition of the key steps of a derivation that highlights 
the close connection between Langevin equations and path-integrals in App. ^ 

To prepare translation into Langevin form we follow [221 and write the RG again in standard form 

^ Z[U] = -V: llvlxly - ilKxly)] Z[U] (23) 



dr 



and define 



K iKxfy = jd'zj^^t.{eum) (24) 

(the "~" indicating adjoint matrices and traces). 

We then abandon the description in terms of weight functionals Zt in favor of one in terms of ensembles of 
fields which are governed by the corresponding Langevin equations. 

Explicitly, to calculate any observable 0\U] of the fields U we write 

{0[U]),^ jD[U]0[U]Z,[U]^l^ 0[U] (25) 

where, separately at each r, the sum is over an ensemble ^[Zt] of N configurations U whose members were 
created randomly according to the distribution Zt. Clearly, for — > cx), the ensemble and Zt contain the 
same information. 

The Langevin equation then schematically'^ reads 

dr [Ua^h = [Ua.it'-h [ Jd'y KlUey]k + <tS] (26) 

where 

is the "square root" of x, x'^y ~ ^xz^zyi and ^ are independent Gaussian random variables with correlators 
determined according to 



>C = |i?K](...)e-'««. (28) 



Note the factor of i which is essential to render this an equation for an infinitesimal change of an element of 
SU{Nc). In fact the components of uj^ := J(Py [S!^y]k[^y]k + can be directly interpreted as the "angles" 
parametrizing a local gauge transformation in transverse space. 

■^This is a continuous time version of the equation that strictly speaking is not unique. The path-integral derivation shown 
in Add. IX] makes it clear that we are to take a "retarded" prescription here in which the derivative on the l.h.s is taken as a 
finite difference and the fields on the r.h.s. are determined at the previous time step. 
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It is of particular importance to note that the possibiUty to formulate the stochastic term via a completely 
decorrelatcd Gaussian noise that is to say with (CS' '^y'"') = S"''^ 6'^^ S^xy , reduces the numerical cost for a 
simulation like this considerably. 

In our practical implementation, wc discrctize the Langevin equation in Eq. H26|) using a regular square 
lattice of iV^ sites (volume (Na)^, where a is the lattice spacing), with periodic boundary conditions. In 
dimensionless units (a = 1) and defining a rescaled, discrete, evolution time 

s := (29) 

the Langevin equation becomes 

U{x; s + Ss) = U{x] s) eyi-p[it''uj°'{x] s)] (30) 

where 

w'^{x,s) = V5'sY,Mx-y)[^-u\x:s)u{y;sT''£!l{y) 

V 

- 5sY^S{x-y)Uy\it'^U\x,s)U{y,s)]. (31) 



Here 



if.(r) = 5, = (32) 



and ^ is a Gaussian noise with {x)^j {y)) = 6ijS°'''Sx,y Ki{r) and S{r) are taken to be = if r = 0. 
In practical numerical work we have to carefully consider the following issues: 

Boundary conditions: The sum over the volume (y) above must be defined in conjunction with the 
boundary conditions used. In general, the periodic boundary conditions minimize the effects of the finite 
volume and makes the system, on average, translationally invariant. Because of the long range of the evolution 
kernel the distance between points x and y must be evaluated taking the periodicity into account. This can 
be done in several non-equivalent ways; we choose one which evaluates the distance {x — y) to be the shortest 
one while taking the periodicity into account. This is achieved by modifying the argument of the functions 
Ki and S as follows: 

{x-y)^ {x~y)+ nN, (33) 

where Ui G { — 1,0, 1} is chosen to minimize the absolute value of the r.h.s. of Eq. H33|l . This procedure cuts 
the range of the evolution kernel to be < N/2, i.e. we take the volume integral over a single periodic copy of 
the volume. 

Let us note here that the boundary condition is significant only if the saturation radius Rg ~ l/Q^ is of order 
system size. As we shall see below, the boundary conditions (or the system size) have a negligible effect on 
the evolution as long as 1/Qs N. 

Fourier acceleration: Due to the non-locality of the evolution operator, one time step evolution of the 
N X N system using Eq. (|31|l requires oc numerical operations. This is impractical except for the smallest 
volumes. However, we note that Eq. H31(l can be decomposed into convolutions as follows: 

cc;" = [c(/f„er) - {U^T'^CiK^, ?7"''e-)j - (Ss) \it'')''^{U^y'^CiS, C/'*^)j . (34) 
The convolution can be evaluated with Fourier transformation: 

C(A, B)x = ^ A(aj - y)Biy) = [J-(A).F(i?)]. (35) 
y 

Thus, using FFTs the cost can be reduced to A^^ log A^, which makes a huge difference. Nevertheless, the cost 
is still large: there arc 192 normal or inverse transformations in each time step (note that Ki and S need to 
be transformed just once). 



11 



typeset with I^TeX on February 1, 2008, 22:08 



Time discretization: The discrete time step {Ss) introduces an error oc ((Ss)^ for eacli update step. We 
cliose Ss so tliat tlic largest values of uj°'{x), which characterizes tlie change in U{x) in a single step, always 
remains smaller than 0.3. This resulted to time steps Ss ^ 5 x 10^^, which is overly conservative; even 
quadrupling Ss had no observable effect on the evolution. 

Spatial continuum limit: It is essential to have control over the finite lattice cutoff effects. In the 
continuum the evolution equations are completely scale invariant; the physical length scale is given by the 
initial conditions. On the lattice we have 2 new scales, lattice size L and -spacing a (which equals 1 in 
dimensionless units above). These scales must be chosen so that the physical scale, saturation scale is 
safely between them; in any case, the effects of the finite IR and UV cutoffs must be checked. In short, it 
turns out that the requirement Rg <^ L is easy to satisfy; a factor of 3-4 is sufficient in order to have a 
negligible effect (i.e. there are no IR problems). On the other hand, wc find that a should be several orders 
of magnitude smaller than Rg for the UV effects to vanish, making the continuum limit very non-trivial. 

Initial conditions: For the initial condition we choose a distribution of [/-matrices with a Gaussian 
correlation function: 

(C/tC/^)(xexp(^^^^^). (36) 

Scale R here is close to the initial saturation scale Rg (depending on the precise definition of Rg). 

For a linear (scalar or vector) field it is straightforward to generate distributions with (almost) arbitrary 
initial 2-point correlation functions, (AqAx) = C^, using Fourier transformations: just let cx \/Ck^k, 
where is white Gaussian noise, P(^) oc exp(— ^^), and Ck is Fourier transform of Cx- Field Ax is obtained 
with an inverse Fourier transformation. 

In our case the degrees of freedom are elements of SU(3) group and this method cannot be applied directly. 
However, we can utilize the approximate linearity of small changes in group elements and build up the 
desired distributions in short steps. The method works as follows: let the desired correlation function be 
{U^JJa) = Cx- Wc start from initial state Ux = I and repeat the operation 

Ux ^ e'^-/^Ux (37) 

n times, generating new suitably distributed adjoint matrices Ax = A'^t'^ for each step. Here A'^ is a Fourier 
transform of A'^ oc {C^Y^^^^, where are independent Gaussian random variables and C^ is the Fourier 
transform of — In Cx (The use of the logarithm is a trick to give the correct asymptotics to the correlation 
function.). 

Thus, the matrices Ux "diffuse" towards the desired state by random walking along the group manifold. The 
number of the steps, n, has to be large enough in order to give us properly randomized {/-matrices. Naturally, 
the method is only approximate, but the deviations from the desired correlation function remain small. 

We also used another, completely unrelated, method for generating the initial state: this time wc start from 
a completely random configuration of [/a;-matrices G SU(3), and make the substitution 

Ux ^ expi^R^W^px « (1 + ^^R^W^rUx , (38) 

for n large enough. is evaluated as lattice finite differences. Thus, the last expression tells us to "smear" 
the matrix Ux with its neighbors n times. For linear fields this procedure yields exactly Gaussian 2-point 
function (UxUby) oc exp[— (a; — y)^/27?^]. This is straightforward to verify by making a Fourier transformation 
of Eq. H38|l . However, because we have to ensure that Ux remains on the group manifold, we have to project U 
back to SU(3) after each step. This is done by finding U' S SU(3) which maximizes ti{WU'), and substituting 
U ^ U' - The projection step makes the initial state generation method again approximate. 

Both of these two methods yield initial ensembles which have almost Gaussian initial 2-point functions. Due 
to performance reasons we used mostly the first method (the latter method requires n 3> i in order to really 
produce long-distance correlations). 
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3.2 Numerical results from JIMWLK simulations 



In our analysis wc use volumes 30'^-512^, with 8-20 independent trajectories (with independent initial condi- 
tions) for each volume. This amount of statistics proved to be sufficient for our purposes; actually trajectory- 
by-trajectory fluctuations in the dipole operators are quite small (i.e. already one trajectory gives a good 
estimate of the final result). 

In Fig. |Slwe show the evolution of the dipole operator, averaged over 8 independent trajectories on a 256^ 
lattice, from interval asT = ashi{l/x) — 0...2.5. Qualitatively, we can observe that the dipole operator 
soon approaches a specific functional scaling form, where the shape is preserved but the length scale is 
shortened under the evolution. This is clearly visible on the logarithmic scale plot on the right: the Gaussian 
initial condition quickly settles towards the scaling solution which evolves by moving towards left while 
approximately preserving the shape. 



256 « =22 



256" =22 




100 



x-y 




Figure 5: Evolution of the dipole operator from the JIMWLK equation on an 256^ lattice, shown on a linear (left) 
and logarithmic (right) distance scale. The general pattern as predicted from the existence of the fixed point and 
drawn qualitatively in Fig. |1| shows up clearly. On the right we see a scaling form emerge: a stable shape of the curve 
that merely shifts to the left unaltered. 

From the dipole correlation function we can now measure the evolution of the saturation length scale Rs- 
There are several inequivalent ways to measure Rs from the dipole function; perhaps the most straightforward 
and robust method is to measure the distance where the dipole function reaches some specified value. This 
is the approach we adopt here; precisely, we define RsiT) to be the distance where the dipole function 
N{{x - yf) = l/iVctr(l - UlUy) reaches a value c: 

N{Rl) := c (39) 

where c should not be much smaller than 1 . The precise definition of course is c-dependent - different choices 
will lead to a rescaling of the value of Rg - and we adopt the convention to set c = 1/2. Such a definition 
can be used both inside and outside the scaling region and will encode some physics information in both. 

In Fig. Elwe show i?s(r), measured from various lattice sizes and initial values for Rs{t = tq). Since the 
initial tq is not a physical observable in the sense that we are not in position to match initial conditions to 
experiment, we are at liberty to shift the r-axis for each of the curves in order to match the curves together. 
In Fig. lathis procedure clearly gives us a section of the R(t) curve, independent of the system size or initial 
size. 
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From RsiT) we obtain the instantaneous evolution speed 



A = - 



d\nR, 
dr ■ 



(40) 



Note that contrary to Rs this quantity is unique and independent of the convention for c in Eq. 139|) inside 
the scahng regime. Outside this uniqueness is lost. 

In Fig.Hwe show A plotted against Rs, measured from several volumes and initial sizes. The curves naturally 
evolve from right to left, from large Rs towards smaller values. As opposed to Fig. [HI in this case the curves 
are fully determined by the measurements. Note that for a scale invariant evolution A should be constant, 
independent of Rs- This is clearly not the case here. We observe the following features: 
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Figure 6: The evolution of the saturation scale Rs 



measured from system sizes SO'^ 
shifted in r to match at large r. 



512 . The curves are 



Figure 7: The evolution speed A — —drlnRs plotted 
against Rs for various system sizes and initial Rs- 



• A starts relatively small but increases rapidly. We interpret this to be due to the settling of the Gaussian 
initial state towards the scaling form. This change in shape is clearly visible in Fig. ^ Such behavior 
typical for initial state effects. This effect should be much smaller if the initial ensemble would more 
closely resemble the scaling solution. In these particular simulations, in which we were mainly interested 
in universal features of the evolution equations at larger r, part of the initial changes also arise from 
the fact that the initial Rs may be too large; i.e. the initial state dipole fimction may be too far away 
from unity at distances L/2, the maximum distance on a lattice. 

• At small Rs (large r), X{Rs) clearly settles on an universal curve, which decreases with decreasing Rs- 
This curve is a function of Rs/ a, i.e. the saturation scale in units of the lattice spacing. Thus, this is 
a UV cutoff effect which should vanish (curve turn horizontal) when a — > {Rs/a — > oo). However, 
we are clearly still far from this case, even on largest volumes. In order to obtain a continuum limit, 
one should perform extrapolation a Q. We do this as follows: first, we note that the maxima of 
\{Rs/a) -curves for each of the runs fall on a universal curve when plotted against a/ Rg at maximum. 
This is shown in Fig. |S1 Linear function does not fit the data well; on the other hand, fitting a power 
law + constant we get a reasonable match (shown in the figure). However, with the steepest part of 
the curve in a region without data it is clear that the extrapolation is not safe. Without a definite 
functional form for the extrapolation our data is not good enough to allow for a reliable continuum 
limit estimate. (Nevertheless, we note that the comparison to BK simulation results, discussed below, 
makes the continuum limit in JIMWLK much clearer.) 
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• On the other hand, the IR cutoff effects appear to be mild: the behavior of X{Rs) curves is largely 
determined by the initial radius Rs{to), independent of the system size L (so far as L is clearly larger 
than Rs). This is clarified in Fig.lHl where evohition of Rs{t) is shown from volumes 30^ - 64^, using 
the same Rs{to) w 6a. Only small initial deviations are seen. 




~ a 2, a t = a log 1/x 



Figure 8: Approach to the continuum limit of A/q 
against a in units of Ra . The values are taken from the 
maxima in Fig. |7| The dashed line on the top shows 
the continuum value obtained from a simulation of the 
BK equation. A continuum extrapolation (continuous 
line) from this alone appears unreliable with most of 
the data away from the steep section. 



Figure 9: The evolution of the saturation scale with 
different IR cutoffs (lattice sizes). Except for a small 
initial difference, the differences are well within statis- 
tical errors. 



Thus, the necessity of having a very large hierarchy of the relevant scales a R, < L prevents us from 
making continuum extrapolations. The cutoff sensitivity problem is also present in the BK equation H19() . 
However, the BK equation is much simpler to study numerically than the full JIMWLK equation: firstly, 
there is only a single degree of freedom, a scalar field instead of a SU(3) matrix field, and secondly the 
evolution equation is fully deterministic, instead of statistical. This enables us to study much larger range of 
volumes and initial conditions with reduced errors. 

Moreover, under certain conditions, the BK equation can be expected to be a fairly good approximation of 
the JIMWLK equation: in Fig. ^| we show the difference of 4-field correlators 

{tv{UiUy)iv{UlU.)) - (tr([/t[/^))(tr([/t [/,)). (41) 

For concreteness, we chose the points so that {x — y) \\ ei, [y — z) |j 62 and |a; — y| = \y — z\. The natural 
magnitude for these correlators is ~ 1. Thus, given the initial conditions we used, we see that the correlators 
cancel to a 1-2 % accuracy, with the violations staying roughly at the same size throughout the r interval 
covered. This seems to indicate that for the type of initial conditions used, the BK equation should give a 
good approximation of the salient features of the evolution. In particular we should be able to refine our 
understanding of the UV extrapolation and augment Fig. |H1 with results from BK simulations in which we 
also use a lattice representation for transverse space to parallel our treatment of the JIMWLK equation. 



3.3 BK and JIMWLK compared 

In order to study the cutoff-dependence and to facilitate the comparison with JIMWLK simulations we study 
the BK equation with two different numerical approaches. 
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Figure 10: Factorization violation for a 4-field correlator. Typical error size is indicated by dashed lines for a curve 
in the middle. They are generically small and indicate that the large Nc limit, i.e. the BK equation will be a good 
approximation as long as the initial condition does not contain strong 1/Nc corrections. 



First, we discretize Eq. H19|) on a regular square lattice with periodic boundary conditions as in JIMWLK. 
Using the translation invariance, we can substitute N^y Nx-y, and set, for example, y = in Eq. H19|) . 
obtaining 

a N f / \ 

drN. - ^ y (^3^ +N.-Nx- Nx-. N^) . (42) 

This equation can be readily simulated on a square lattice. The evolution kernel again can be written in 
terms of convolutions, and can be efficiently evaluated using FFTs. We again make the choice that when 
z ^ or z = X, the integral kernel evaluates to zero. The periodicity is taken into account analogously to 
the JIMWLK equation. 

The second method we use has been developed by [S^, similar results had earlier been obtained in [30]. By 
first Fourier transforming the evolution equation (|19|l to momentum space, it can be further transformed into 
a form which depends only on logarithmic momentum variable ^ ~ In fc^. In this form the kernel involves only 
1-dimensional integral, and it becomes easy to include a very large range of momenta in the calculation; in 
our case, more than 20 orders of magnitude. Thus, this enables us to obtain a cutoff-independent continuum 
limit result for the evolution exponent A. The key result we obtained with this method was already shown 
in Fig.|H|for comparison with the JIMWLK extrapolation. 

In Fig. ^2 we show the evolution of N measured from 4096^ lattice. The similarity with the dipole function 
of JIMWLK equation is striking, see Fig.|31 Again we see the initial settling towards the scaling form, which 
evolves by shifting towards smaller R while preserving the shape. Also here we see that, even with a scaling 
form reached, the spacing of curves and hence A is not constant. This again is due to the growing influence 
of the UV cutoff. 

We are now in a position to complement our results from the JIMWLK equation. We extract A in the same 
way as in Figs. [T] and |S1 This time, there are no statistical errors in the measurements of A. A function of 
form A/tts = ci +C2{a/RsY is an excellent fit to the data, with result ci = 2.26, C2 = 2.0 and = 0.32. The 
a ^ result agrees perfectly with the momentum space result. This is shown in Fig.^] The agreement of 
BK and JIMWLK results at corresponding cutoff scales is remarkable over the whole range. Most striking 
is an evaluation of the active phase space. While on the IR side activity is limited to within one order of 
magnitude of the saturation scale Qs, on the UV side we need the cutoff about 6 orders of magnitude larger 
than Qs to get within 1% of the continuum value. as is easily extracted from the power law fit. 
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Fixed coupling, volume 4096 





r I a 

Figure 11: The evolution of N{r) on a 4096^ lattice, 
measured along on-axis direction. The curves are plot- 
ted with interval 5t = 0.25/a, and the total evolution 
time is 7/a. 



Figure 12: Continuum extrapolations for A using 
JIMWLK and BK simulations. The UV cutoff is plot- 
ted in units of the physical scale Rs. The continuous 
line is A/a^ = 2.26 - 2.0(a/i?s)°'^^ 



This situation is clearly unnatural: large open phase space of this magnitude usually generates large loga- 
rithmic corrections that need to be resummed to get physically reliable results. In the following we will argue 
that the main source for such corrections are due to running coupling and study the effects which, among 
other new qualitative features, will bring a large reduction of active phase space in the context of the BK 
equation that erases the potential for large logarithmic corrections. 

4 Running coupling effects 

To account for running coupling effects in the simplest possible way, wc decide to have the coupling run at 
the natural scale present in the BK-cquation, the size of the parent dipole {x — y)^, i.e. in Eq. (|19|l . we 
replace the constant as by the one loop running coupling at the scale l/{x — y)^: 

a,^a,{l/{x-yf) = — ^ with /3o = (lliVc - 2iV/)/3 (43) 

A in Eq. (|43|l is related to Aqcd by a simple factor of the order of 2tt. Such a simple replacement, of 
course, is not what we expect to arise from an exact calculation in which the two loop contributions that 
generate the running will emerge underneath the integral over z. This comes on top of the fact that at next 
to leading order there will be corrections beyond the running coupling effects that are entirely unaccounted 
for by this substitution. Nevertheless, we would like to argue that Eq. H43|) . when used in the BK equation 
captures the leading effect associated with the "extraneous" phase space uncovered in Sec. 13.31 This has 
been the outcome of a long discussion in the context of the BFKL equation and its next to leading order 
corrections |34ll35irT()ll37ll38ir^l4()| Here we simply take this as our starting point with the attitude that a 
scale setting technique like the one employed in Eq. H43(l is but the easiest way to get even quantitatively 
reliable initial results. We plan on supplementing this with more detailed studies based on dispersive methods 
of implementing the running coupling effects. 

With running coupling wc have to face the presence of the Landau pole in Eq. (|43|) . The main point here is, 
that while indeed we do not control what is happening at that scale wc are justified in providing a prescription 
for its treatment. With Qs ^ Aqcd the integral on the r.h.s. of the BK equation yields extremely small 
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contributions at scales of the order of Aqcd, typically of the order of 10~^ times what we see around Qs- 
This is the reason why we are allowed to ignore whatever happens there and to freeze any correlators such 
as N at these scales at the values provided by the initial condition - a reasoning that is already necessary 
at fixed coupling in some sense. If we turn on running coupling then, the physical thing to do is to switch 
off the influence of the divergence at the Landau pole, for instance by freezing the coupling below some scale 
/iQ. Of course we make sure that the details of where we do that will not affect our results. 

Such an argument about the integrals on the r.h.s. is valid at 5 « or infinite transverse target size where the 
nonlinearity is present and sizable everywhere we look. For large h one necessarily enters the region where 
the nonlinearity is small and there one would remain sensitive to any value of /xq. That is our reason for not 
considering such situations. 

Let us now explore how this modification affects evolution. The crucial first observation is a clear reduction 
of active phase space as shown in Fig. 1131 on the left panel we show an example for the extrapolation of 
A/q:s(Qs(t)^) to its continuum value at one value of t (corresponding to as{Qs{T)'^) = 0.2) and compare to 
the fixed coupling case. With running coupling we are able to reach the continuum limit already for relatively 
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Figure 13: Left: example of a continuum extrapolation of \/aa{Qs{TY) for a fixed initial condition at a r-value 
where Qs(Qs(r)^) = 0.2, compared to the fixed coupling case. Right: A(i?s(r)), measured using different lattice 
cutoffs, together with the extrapolated continuum limit curve. The dashed vertical line indicates the point where the 
values on the left are extracted. 



small lattices, i.e. with a UV cutoff that would be orders of magnitude too small in the corresponding fixed 
coupling run. Note that as a consequence the continuum value is reduced by a noticeable factor. On the right 
panel we demonstrate how the extrapolation was done: for a given initial condition we perform runs with 
different lattice spacings and measure A(t) and Rs{t). We get a reliable continuum limit extrapolation for A 
in a wide range of r-values; when Rs becomes small 5-lOa) the measurements start to fan out. However, 
as will be shown in Fig. 1151 we can extend the measurements by starting from different initial i?s-values. 

We will turn to a more systematic study below and only emphasize that the results given in what follows are 
physical in the sense that they apply to the continuum limit. 

First and foremost, the running coupling reflects the breaking of scale invariancc. It would be natural 
to expect that the dipole function now might never exhibit the scaling behavior of Eq. H21a|) that was so 
characteristic at large r for fixed coupling. In fact numerical evidence suggests otherwise. As shown in 
Fig. I14al which shows results Nr^xy for large r against rQs{T). The different curves only slightly deviate 
from each other: scaling is only violated very weakly. We attribute this to the fact that Qs{t) ^ Aqcd and 
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'■/R/x) a T = a In l/;c 

Figure 14a: Approximate scaling behavior of N^^^y Figure 14b: Scheme and t dependence of A for one 
on the asymptotic line of Fig. [HI initial condition 

thus the scale breaking is not strongly communicated to the active modes. Where the form of Nr.xy shows 
approximate scaling it makes sense to talk about Qs in complete analogy to the fixed coupling case. Again 
defining the rate of change A according to 

A(t) := 5, In Q,(t) , (44) 

leads to an almost unique and scheme independent quantity. The main new feature is that now, in contrast 
to the fixed coupling case, A is r-dependent. Both scheme independence in the near scaling region and r- 
dependence of A arc shown for one initial condition in Fig. I14bl The curves correspond to different values of 
c in Eq. (|39|l and a direct definition of A via a moment of the evolution equation to be given below. We see 
first a strong rise which occurs in the pre scaling region that turns into a slow decrease in the "near" scaling 
region where the curves converge to a single line. The accuracy to which the the curves agree in the near 
scaling region at large r can be used as a measure of the quality of the scaling. 

Since in the scaling region, Qs{t) is the only scale left in the problem, it is natural to plot A(t) against 
Rs ^ 1/Qs- This allows to compare the asymptotic behavior of simulations with different initial conditions 
as shown in Fig. El In this plot we find that all A values obtained from Gaussian initial conditions are 
limited from above by an asymptotic curve that falls with growing r (from left to right in the plot). It is on 
this asymptotic curve that the scaling shown in Fig. I14a| holds. 

With the global features established let us turn back to a quantitative analysis of active phase space and a 
study of the t dependence in the near scaling region. 

To this end we first show a detailed analysis of the continuum extrapolation of one trajectory at several 
values in r in Fig. 1161 On physics grounds we expect the UV cutoff needed in a given case to be determined 
by Qs at that the r considered. In Fig.^l left, we sec that with growing r the extrapolations of A towards 
the continuum limit at the left reach the asymptotic plateau later and later. That this indeed follows Qs is 
shown in the middle. What remains open in these plots is the precise form of the r-dependence of A. The 
numerical result here is surprising: our simulations show a clear scaling of A/ as{Qs{T)) in the asymptotic 
region over the range explored in this simulation. This is shown in Fig. 1161 right, which summarizes our 
results on t and scale dependence for A in the range of Qs values appearing in this simulation. 

Active phase space is clearly centered around Qs{t) both in the UV and the IR. The reasons for the bound- 
edness are very different on the two sides: while the IR safety is induced by the nonlinear effects already 
present at fixed coupling, the limited range into the UV is solely caused by the running of the coupling. 

While it is natural to express the r dependence of A via Qs (r) , it is rather surprising that for the small values 
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Figure 15: A = dT-lnQsir) for different initial conditions. All curves starting from Gaussian initial conditions stay 
below an asymptotic line along which we find (approximately) scaling dipole functions A'^. 
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Figure 16: UV dependence of the BK-equation on the asymptotic line after the parent dipole scale is taken to 
determine the coupling. Left: A against a in units of lattice size. Middle: A against a in units of R^, the natural 
units. Right: the same rescaled by YQs(Qi). Active phase space is limited in the UV (and IR) by 5 to 10 x Qs. The 
r-dependence of A is, to a good approximation, given by A(r) — Xq^ asiQsir)'^). 
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of the coupling constant in the range covered in Fig. 1161 that an expansion of A of the form 

A(t) ~ casiQsi'T')'^) + small corrections (45) 

appears to be totally inadequate and instead a fit of the form 

A(t) = Xo^asiQsir)^) (46) 

appears to be surprisingly successful. If the latter provides a good fit, then the corrections in Eq. (|45|l must 
be important in the range considered. As it turns out, the fact that A(r) can not be expanded in a simple 
power series in as{Qs{T)^) is linked with the large reduction in phase space compared to the fixed coupling 
case. To see this, let us use an idea by McLerran, lancu, and Itakura |3T] that was introduced in the fixed 
coupling case to show that there A is necessarily constant in the scaling region. 

They assume (in the context of the fixed coupling BK equation) that Nj- depends only on the combination 
(a; — y)^Qs{T)^ =■ i"'^Qs{t)'^- Then r-derivatives can be traded for derivatives and the l.h.s. of the BK 
equation is easily rewritten as 

drN{r^Q,{T)^) = r^dr2N{r^Qs{T)^)drliiQs{T)^ = r^dr2N{r^Q,{T)^)2Xr . (47) 

Taking into account the spatial boundary conditions imposed by "color transparency", the fact that N{r'^ = 
0) = and N{r'^ = oo) = 1, leads us to take the zeroth moment of the BK equation -to integrate with 

.2 

in order to isolate A on the l.h.s.. The r.h.s then provides an integral expression for it in form of its 
zeroth moment: 

2-A(r) = ^ / ^^iN{n'Ql)+N{v'Ql)^N{v'Ql)^N{n'Ql)N{v'Ql)) (48) 

where we have set u = x — z v ^ z — y ioi- compactness. Due to scale invariance of the integral, the r.h.s. 
is indeed independent of r: the scale common to all N can be chosen at will without changing the integral. 
Indeed, in the scaling region Eq. H48|l is completely equivalent to Eq. H44|l . Outside that region it can be used 
as an alternative definition of a quantity that gives a qualitative description of the rate of evolution with the 
same justification as any other definition of A or Q^. Its main advantage is that now it is directly related 
to the integral expression driving the evolution step and we can directly analyze active phase space in this 
expression where needed. 

On the asymptotic line with near scaling the argument is only approximate, but to good accuracy we obtain 
2-A(r) = §| ^^a4iy){N(u'Ql) + Niv'Ql)-Niv'Ql)-N{u'Ql)N{^^^^^^ (49) 

and we can take the r.h.s. as a definition of A(r), ignoring the small corrections. Fig. I14b1 shows that we 
indeed agree with the conventionally measured quantity. 

Due to the near scaling of iV also in the running coupling case together with all other ingredients of the integral, 
it is really only the presence of the running as underneath the integral that distinguishes the integrals in 
Eq. H48|l and Eq. H49() . To obtain the strong relative reduction in active phase space observed in Fig. E| we 
need as to vary considerably over the range of scales contributing to the integral in the fixed coupling case. 
Wherever this happens, A(t) will be a nontrivial function of as{Qs{T)^)- Conversely, at asymptotically large 
Qs, where the as(l/r^) probed in Eq. H49|) turn essentially constant over a range comparable to the phase 
space needed in the fixed couphng case, we should return to a situation as in Eq. The price to pay is 

the reopening of phase space to the size encountered already at fixed coupling. 

To illustrate this we list active phase space in terms of a/i? needed to find A within a given percentage of its 
continuum value at different as(Qs(T)^) in comparison with the fixed coupling result. 

In the far asymptotic region we arc back to the unsatisfactory situation that we most likely miss important 
corrections with the additional drawback that we now have no real idea of what those might be. Fortunately, 
this would appear to occur only at such extremely large values of Qs for this to be purely academic. 
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fixed coupling 



a/i? = 0.0011 
a/i? = 0.00013 
a/i? = 8.2-10-^ 

Table 1: Phase space needed to extract continuum values for A is reduced for larger values of the coupling. The fixed 
couphng situation reemerges asymptotically. 
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With this caveat, we are able to understand the behavior of A on most of the asymptotic line. Fig. 1 171 shows 
the behavior of A as a function of as(Qs(T)^) and Qs demonstrating in the center the large range over which 
a fit according to Eq. H46|) is successful. To the left of this, we see a deviation that would appear to be 
consistent with a turnover towards a linear behavior in a^ . Indeed a purely phenomenological two parameter 
ansatz of the form 

2.26a,(Q,(T)2) 



A(r) 



(50) 



l + aa,{Qs{TYf 

provides a fit consistent with all the available information. Here 2.26 is the value obtained from the fixed 
coupling BK-simulation which must emerge due our argument above and a and h have been fitted to be 6 
and 0.8 respectively. 
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Figure 17: A(r) as a function of as{Qs{'ry')- The central region of the asymptotic line where the running of the 
coupling reduces phase space is best fitted by a A(r) = Aoya7(Q7(7p) or equivalently by the use of Eq. 15611 . At 
large Qs the curve turns towards a form that allows a perturbative expansion of the form A(r) = CQa(Qs(T)^) + 
small corrections at the price of a reopening of phase space. 



To establish the link to the asymptotic behavior given in ^2E21j we first observe that the small coupling 
limit of Eq. (|50|l with its leading 



A(t) = 2.26as(Qs(''')^) + small corrections 
is in full agreement with the asymptotic behavior derived in j41ll42| . This we render here as 

.90 0.47 



A = 



VFTT^ (r + ro)5/6 



higher inverse powers. 



(51) 



(52) 
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where Yq stands for constants determined by initial conditions.^ To show agreement with the first term in 
Eq. 1)52(1 analytically, let us consider a generic power like dependence of A on as{Qs{T)^) to incorporate both 
Eq. ((51|l valid at asymptotically small as and Eq. (|46|l which gives a simple representation for the bulk of 
scales. We make the ansatz 



Xr := dr HQ s{t)/ Aqcd) - Aq [a,(g,(r)2)]' 



47r 



/3oln(g,(T)2/AQCD) 



with l3o = (llA^c - 27V/)/3, and find 



^^QCD 



cxp 



{[< 



Recovering A from this we find 



A(r) 



l)2Ao(|^)"(r-ro) 



\ / 47r 



[ln(- 



^^QCD 



♦ l,Ao^2.26 



(53) 



(54) 



(n + l)2Aoyf(r-ro)+[ln(%g^)] 



1 n+l 



(55) 



so that the leading term of Eq. 1(52(1 corresponds to rt = 1 as expected. 

To extend the comparison beyond that wc simply integrate A in Eq. 1(52(1 back to get an expression for Qs 
and fit the integration constant c to match our simulations. [Alternatively we could have directly used the 
expressions for Qs in (41[I42( and adjusted the unknown parameter.] We find that 



A 



QCD 



c • exp[2 • 0.9(r + Yq)^ - 6 • 0.47(t + Yq 



(56) 



with c ~ 1.5 leads to excellent agreement over the whole range. This is shown in Fig. 1171 

To summarize the results of this section we state that Fig. 1161 also shows impressively that the contributions 
to A arise from within little more than an order of magnitude of the saturation scale [the region in cutoffs over 
which there is an appreciable change in the value extracted]. This indicates that the dominant phase space 
corrections are indeed taken into account and the basic physics ideas that had initially motivated the density 
resummations are realized. While the IR side of the resummation was working already at fixed coupling, the 
UV side of phase space is only tamed if the coupling runs. Only with both the nonlinearities and a running 
coupling are we justified to claim that the evolution is dominated by distances of the order of Rs and thus 
the coupling involved in typical radiation events during evolution is of order ctsiQsi'r)'^). Fig. llTl shows that 
on the asymptotic line A(t) = Ao y^3?7(Qs(Tp) over most of the region where running coupling limits phase 
space. 



5 Evolution and 74-dependence 

It has been argued in the context of the McLerran-Venugopalan model that going to larger nuclear targets 
(larger nuclear number A) ought to have an enhancing effect on the importance of the nonlinearities encoun- 
tered. If we are at small enough x (high enough energies) for the projectile to punch straight through the 



target going to a larger target means interaction with more scattering centers as indicated in Fig. 2(b) Tak- 
ing these to be color charges located in individual nucleons inside the target nucleus the interaction will occur 
with more and more uncorrelated color charges the larger the target. Viewed from up front this amounts to 
a smaller correlation length in the transverse direction. In terms of Qs this would lead to 

Qf(r)2=^^Q^(r)2 (57) 

without a clear delineation on the r values for which this would be appropriate. This simple argument 
has been uniformly used in 0^1 to relate cross section- (or F2-) measurements on protons at HERA with 

^This form may be obtained from a Y derivative of Eq. (83) in 1411 or from Eq. (53) of 1421 . We have converted their Y to 
our T and corrected for a factor of 2 in the definition of A 
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measurements on nuclei at NMC and E665. Despite the fact that the nuclear data used are at rather small 
and large x the agreement found was surprisingly good. 

In the context of small x evolution such a statement becomes less trivial: the whole idea of deriving a small 
X evolution equation is based on premise that the logarithmic corrections determining the evolution step 
are entirely target independent. The only place where target dependence is allowed to enter is the initial 
condition. 

In simple cases, these two notions are not conflicting with each other: indeed a rescaling like in Eq. (|57|1 is 
fully compatible with fixed coupling evolution on the asymptotic line, where (in the continuum limit) 

Q,(t) = e^^Q,(To) (58) 

with a universal value for A, so that a simple rescaling of the properties of the initial condition according to 

(ro)2 = A^Ql(T^f (59) 

automatically yields Eq. (|57|l at all r declaring it fully compatible with evolution. 

At running coupling, however this ceases to be the case. The straightforward physical argument is that 
evolution is characterized by as((5^) which initially for protons is larger than for nuclei because of the 
difference in saturation scales at the initial condition. In this situation evolution will be faster for smaller 
targets and eventually catch up with that of the larger ones. 

Technically we find that with running coupling (again restricting ourselves to the asymptotic line), we find a 
r dependent A that we may parametrize through some function / of ois{Qs{j^'^^ as as in the previous section: 

A(t) - a,Q,(T) = /(«,s(Q.(t)2)) (60) 

In this generic case 

Q^(^)^g/;'^-7(..(Q.M^)Q^^^) (61) 

and a rescaling of Qf (tq)^ will not simply factor out as in Eq. H58() . Therefore and a simple rescaling like 
Eq. (|57|l will conflict with evolution. All phenomenological fit functions encountered in the previous section 
share this property which is at the heart of the reasoning used by Mueller in 03] to argue that evolution will 
eventually erase A-dependence in Qs- 

To illustrate the point, let trace the the memory loss with a simple parametrization like Eq. H53|l . equivalent 
to a T dependence of Qs according to Eq. 1)54(1 . 

Scaling in the A-depcndence at tq according to Q^(to)^ ~ A^Qf (tq)^ within Eq. H54|) then predicts both r 
and A dependence of the saturation scale. 

We find that evolution will indeed strive to wash out A-dependence in qualitative agreement with 01]. Indeed, 
the ratio of saturation scales in the bulk of the region shown in Fig. El then behaves like 



exp 



exp 



{n + l)2Ao(f )"(r - r„) + [ln(^ig£^)]"+^ 



(n + l)2Ao(f)"(r-ro)+[ln(%g^)]"+^ 



1 . (62) 



While the parametrization with n = 1/2 is not applicable at arbitrary large r, this formula nevertheless gives 
a reasonable estimate of the rate of convergence within its window of validity. This is plotted in Fig. 1181 
Clearly larger nuclei suffer stronger "erasing." It should be kept in mind that realistically one should not 
expect real world experiments to cover more than a few orders of magnitude in x and hence too large an 
interval in r — Tq. 

The fact that Eq. (|46|l and the ensuing r dependence is only valid on the asymptotic line of Fig. 1151 prevents 
us from making a direct statement about current and future experiments. Such stronger statements require 
a more thorough study of initial conditions in direct comparison with existing experiments. If dependence on 
initial conditions can be clarified A dependence can be studied numerically also in the pre-asymptotic region. 



24 



typeset with l5rEX on February 1, 2008, 22:08 



Q /GeV 



6 



10 



100 




1 



proton 











5 



10 



15 



20 



25 



30 



X - 



'0 



Figure 18: ^-dependence for proton saturation scales from 1 to 100 GeV. 



6 Conclusions 

In this paper we have, for the first time, shown results from a direct numerical implementation of the 
JIMWLK equations. The simulations have been performed using a Langevin formulation equivalent to the 
original equations. These then have been implemented numerically using a lattice discretization of transverse 
space. 

Our simulations have been mainly concerned with extracting universal features and asymptotic behavior in- 
duced by the equation itself without direct comparisons with data and systematic studies of initial conditions. 

We have clearly established infrared stability of the solutions: once provided with an initial condition that 
furnishes a (density induced) saturation scale, or correlation length Rs ^ ^/Qs, the solutions are completely 
insensitive to IR cutoffs sufficiently below that scale. Such behavior had already been found in the large Nc 
limit in the context of the BK equation and here merely serves as a consistency check. 

Simulations of the BK equation have up to now been only carried out in momentum space, our simulations 
are the first in coordinate space and thus we are the first to show the approach towards the fixed point at 
infinite energy with vanishing correlation length directly. 

Beyond confirming analytical results already found in |22| . we also see the emergence of scaling solutions, 
extending earlier results from simulations of the BK equation to finite Nc- Within the class of initial conditions 
studied these emerge universally already at finite energy. 

The initial conditions studied were all transversally homogeneous and did not contain large intrinsic Nc 
violations. We found that JIMWLK evolution did not enhance these to above the percent level. In such a 
situation the BK equation should provide a good approximation. One should be cautious to generalize this 
too far: any type of lumpiness in the initial condition, hot spots of any kind, would presumably spoil this 
feature and the factorization limit would not be adequate. 

Simulations of the BK equation had already clearly shown that evolution populates a large range in phase 
space of tens of orders of magnitude, a fact easily attributable to the random walk property of the BFKL 
equation which due to density effect is tamed in the IR, but not in the UV. While previous studies have 
taken this as a given, we were forced by the necessity to work on a coordinate lattice to carefully study the 
UV properties to obtain a valid continuum extrapolation. 
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Already within the JIMWLK simulation we again found evidence for enormous active phase space that 
induced relatively large errors in the continuum extrapolations. Since our initial condition did lead to 
simulations which respect factorization we were in a position to refine this study by implementing the BK- 
equation in coordinate space with the same type of lattice regulator used in the JIMWLK simulations. We 
found excellent agreement between the JIMWLK and BK simulations on finite lattices and were able to 
confirm that the (constant) evolution rate of the saturation scale on the asymptotic line A ~ drQs{T) gets 
contributions from the far UV, up to 6 orders of magnitude away. In such a situation one typically expects 
large logarithmic corrections that need to be resummed. Indeed, from the experience with the BFKL equation 
at two loop, there was an obvious candidate for the main effect, the running of the coupling. 

This leads us to the second part of the paper. Here we studied the effect of a naive implementation of 
running coupling on the same type of questions as studied above. This has been done in the context of the 
BK-equations for reasons both of efficiency and principle. The latter will be briefly discussed below. 

The main effect of running coupling was to drastically improve the convergence towards the continuum limit 
over a gigantic range of Qs values between 1 and well above lOOGeV. [Note that the Golec-Biernat+Wiisthoff 
fits to HERA data have rise from 1 to 2GeV from x = 10"^ to lO"'^.] 

Despite the explicit breaking of scale invariance we find that the dipole function still converges towards a 
scaling form with all the t dependence taken by the saturation scale Qs{t)- 

The evolution rate is effected more strongly: A(t) = drhiQair) now is necessarily r dependent. Within 
the window of strongly reduced phase space between 1 and lOOGeV it is well described by a simple fit of 
the form A(t) = .37 as{Q s{t)^) . This is strongly different from the naive expectation that one should find 
A(r) = cas(Qs(T)^) + small corrections. We have numerical and analytical indications that such a behavior 
would reemerge at extremely large r at the price of a reopening of phase space. The explanation for this is 
simple: as shown in Sec. 01 the running coupling appears in an integral expression for A that reduces to the 
corresponding fixed coupling expression where the running is weak over phase space active in the integral. 
Thus a strong reduction in phase space necessarily means that can not be factored out of the integral. At 
extremely large Qs, however, the running will slow down sufficiently to allow the fixed coupling features to 
reemerge. The agreement with |41II42| on the other hand is excellent. 

The running coupling induced t dependence of A is responsible for another feature of small x evolution 
already discussed by Mueller the washing out of the A dependence of Qs inherent to the McLerran- 
Venugopalan model. Scaling Qs at the initial condition at small tq according to the McLerran-Venugopalan 
rule Qf{To)'^ = A'sQp{to)'^ before evolution adds gluons in an obviously A independent manner leads to 
a fully determined r and A behavior of Qs in which, because of the smaller as at larger scales evolution 
of smaller targets can catch up with that of larger ones. Our simulations give quantitative results on the 
asymptotic line. Quantitative studies with specific initial conditions are possible. 

The studies presented in this paper have shown the feasibility of JIMWLK simulations but at the same 
time have highlighted the necessity of implementing running coupling corrections to obtain reliable results. 
In particular the latter has been done with a naive scale setting for the running coupling at the parent 
dipole size. Here clearly improvements are necessary. A more refined treatment would most certainly lead 
to modifications of our quantitative results such as the fits obtained for the r-dependence of A, although we 
expect our qualitative conclusions to remain valid. An example that clearly demonstrates the insufficiently 
of the scale setting used is the fact that such a treatment would not be possible is that such a scale setting 
is impossible in the full JIMWLK context. There the closest thing to a parent dipole running would be to 
set the scale at l/{x — y)^ in Eq. 0. This, however, is incompatible with BK parent dipole running as is 
would simply set the scale to infinity in the a terms of Eq. H23|l - a step incompatible with the derivation 
of the BK-equation. Preliminary studies together with Einan Gardi and Andreas Freund have shown that 
consistency in this respect can be achieved if one implements running coupling with dispersive methods (as 
widely used for other reasons in the renormalon context), but further work is necessary to gain thorough 
understanding of all the issues involved. 

Complementary issues like the proper choice of initial conditions and questions of 6-dependence need to be 
studied to be able to put comparison with experiment on a sound basis. 
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A From Fokker-Planck to Langevin formulations, the principle 

Although this topic is covered in many textbooks, wc find that the presentation is often somewhat oldfash- 
ioned and the intimate connection with path-integrals is not always mentioned. Let us therefore begin by 
demonstrating the relationship of Fokker-Planck and Langevin formulations with the aid of a simple toy ex- 
ample, a particle diffusion problem, that, for ease of comparison, shares some of the features of our problem. 
The toy model equation is 

drPr{x) = -HfP Pr{x) (63) 

with a Fokker-Planck Hamiltonian 

HFP^lidlixMx)^d: . (64) 

The particles coordinate x in D dimensional space corresponds to the field variables of the original problem 
and the probability distribution to Zr- It serves to define correlation hmction of some operator 0{x) via 

{0{x))r ^ Jd^xO{x)Prix) 

Eq. H63() admits a path-integral solution which is, as usual, built up from infinitesimal steps in r: this is a 
solution P{t' ,x';t,x) to Eq. H67|l. with initial condition 

P{T,x';T,x)^d'^°\x' -x) (65) 

that is composed of many, infinitcsimally small steps at times t + fee. In this solution, the product rule 



P{t',x';t,x) ^ J d^Xn^i . . .d°xi P{t' ,x';Tn-i,x)P{Tn-i,Xn-i;Tn-2,Xn^2) ■ ■ ■ P{ti,Xi;t,x) (66) 

expresses the finite t' — t solution P{t' , x'\ r, x) in terms of infinitesimal steps P(T.i, Xi] r^^i, aJi-i). Solution 
with arbitrary initial conditions Pro{x) are then recovered via Pt{x) = / (i^yP(T, a;; tq, y)PTo (?/)• 

The derivation of the infinitesimal step from r^-i to = Ti-i -I- e is textbook material and wc omit the step. 
It reads (the index i on the coordinate also refers to the time step, not to the vector component) 

Pin,X.,;T.,^l,X,^,)=N jd^p, g-e(ipfpfx°''(-.-l)+^<P.-(2^^~-(x.-i))) +0(,2) (g7) 

where <y^{x) ^id^Xafi{x) and N a coordinate independent normalization factor. 

From here, the step to the Langevin formulation is trivial: as this expression is quadratic in the momenta we 
can trivially rewrite it with the help of an auxiliary variable [i again the step label] as 

P{n,x,-n^i,x,^^) ^ N jd^'i, Vd^tW^e-^«''^"'(==-)«- 6''{x,-[x,^i+e{i, + a{x,^i)]) . (68) 

The (5^(. . .) arises from the momentum integration and determines Xi in terms of Xi^i and the correlated 
noise ^^.^ The equation for Xi, 

X, = Xi^i + e{ii + (7{xi^i) , (69) 

is called the Langevin equation. It contains both a deterministic term [the a(xi_i) term] and a stochastic 
term [the ^ term]. To fully define the problem without any reference to its path- integral nature one then 

^Rc-exprcssing the delta function via a momentum integral and performing the Gaussian integral over 4i immediately recovers 
Eq. 
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needs to state the Gaussian nature of the (correlated) noise separately. This is often done in textbooks. We 
find that the path-integral version is much more elucidating.^ 

If, as in our case, x factorizes as 

Xfivix) = £i_ia{x)£i,a{x) , (70) 

a simple redefinition of the noise variable leads to a version of the Langevin formulation with a Gaussian 
white noise: 

P(T„a;,;;T,;_i,a;,_i) = N Jd"(, e'^^' S''{x, - [a;,_i + e(£:(a;,_i)e, + a(a;,_i)]) . (71) 

The correlation is now absorbed into the Langevin equation^ which reads 

X, = x,_i + e{E{x^^i)i, + ct{x,_i) , (72) 

Wherever this form is available, it will be the most efficient version to use in a numerical simulation as the 
associated white noise can be generated more efficiently then the correlated noise of the more general case. 

What is left is to sketch a numerical procedure to implement the solution given in Eqns. (|66|l . (|67|) using 
the expressions Eq. (|68|l or Eq. H71fl . As a first step we replace the average with the probability distribution 
Pr{x) by an ensemble average 

{0{x))r ^ I d^x 0{x) Prix) w ^ 0{x) (73) 

ueE[p^] 

where the ensemble is taken according to the probability distribution Pt{x). This we can easily do for the 
initial condition at tq. The Langevin equation then allows us to propagate each ensemble member down in 
r by e, thereby giving us a new ensemble at ti = r + e. Iteration then gives us a chain of ensembles for 
a discrete set of that provide an approximation to a set of P,-. (a;) in that they allow us to measure any 
correlator {0{x))r according to Eq. (|7^ . 

The above clearly shows the nature of "Langevin system" as a way of writing a path-integral solution to a 
differential equation that is particularly easy to implement. Many often confusing issues, such as the time 
step discretization issues often discussed, find their natural resolution in this setting. 
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